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ABSTRACT 


Four time differencing schemes were tested using a 
barotropic primitive equation model on a spherical 
staggered grid with an analytic input in order to com- 
pare amplitudes, phase speeds, and computation time for 
each. The methods tested were the leapfrog, Euler- 
backward, leapfrog-trapezoidal,» and Adams-Bashford. One 
set of experiments was performed uSing an averaging 
technique to reduce the effects of gravity waves in the 
higher latitudes. Another set was performed without the 
averaging in order to determine the effects of this 


technique on the solutions. 
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ie ee hODUCT ION 


In the field of operational numerical weather prediction, 
the trend, in recent years, has been toward the diese ila oneine 
Seesophnisticated global prediction models. This has been 
made possible by the rapid expansion of computing capacity 
and developments related to general circulation research. 

The purpose of this study was to examine various time 
differencing methods, using a barotropic primitive equations 
model on a global staggered grid. A spherical harmonic 
analytic stream function was used for the initial conditions. 
By using an analytic initial condition, errors in real data 
observations and analysis, which are unavoidable in practical 
dynamical prediction, were eliminated. 

The objective was to compare the time required for com- 
putation, amplitudes, and phase speeds for each of the time 


differencing schemes. 
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Ii. 


BaAROURor LOePREMELTEVE EQUATIONS MODEL 


Time differencing experiments were performed using the 


free surface barotropic primitive equations. The integra- 


tions were carried out on the sphere using the difference 


method of the Arakawa type which was developed by Winninghoff 


meo7i). 


mo FRIMITIVE EQUATIONS 
The primitive equations, in 


flux form, for this model are: 








3 (uh) i 3d (uuh) 
= - ———— |— ++ 
at a cos 0 ]|dA 
uvh tan 0 
ee + fy = 
a 
d (vh) 1 d (vuh) 
== | — + 
dt a cos O] dA 
uuh tan 0O 
- + fuh - 
a 
d (h) 1 0 (uh) 
et a cos O]| 0A 


spherical coordinates and in 


d({uvh cos 6) 


00 


h (gh) 
(1) 





a cos © dA 


d(vvh cos 0) 


d © 


h d(gh) 


a or 


(2) 





d(vh cos 606) 
+ 
a0 


(3) 


Equations (1) and (2) are, respectively, the zonal and merid- 


ional momentum equations, and equation (3) is the continuity 


equation. 


el: 





B. GRID 

The spatial finite differencing was performed on a 
staggered, spherical grid. The wind and height variables 
were carried at alternate points (see niga 1) with height 
only at the poles. The latitudinal and longitudinal grid 
increments were five degrees. This gives 2560 points 
(72 x 35) over the globe, with wind and height carried at 


1260 points each. 


h (North Pole) 


Mey ht Ugve oh SEs: = as: |OOR 
h Vv. h Vs cee rs « « U,V 
ose h Es h Ce ceote: csc ce 
( h gan h RVs omens. s20s « U,V 
=] vas lat Bl Ni h Mee cots. « 
ia h (South Pole) 


FIGURE 1. LOCATION OF VARIABLES 


OF SPATIAL PONTE DEE Ee RE NCING 

Spatial differencing of the Arakawa (1966) type was 
used which eliminates the spurious energy growth which can 
occur with standard finite difference approximations to 
the nonlinear advection terms. The difference equations 


used here for this model are: 


me, 





Cen k 
A (uh) 33 1 eee cuneate at) (ah) j 1, 














At a cos 0. OAX 
a (u; 4+Uj 442) (vh) i541 00S Sh EE oN 3) (vh) | 4.1. ¢0s O, 4 
2A0 
u.s(vh):. tan 0: Bee ae lh = hh... 
+ aos tT ae £. (vh) ; _ eee se 14 me kg (4) 
s td a cos 05 ne 
A(vh) 35 7 I ees pee cee ye) (uh) 5 on 
At ancOs O- es hee i — 
§ (vs 5*¥i5e+2) (vb) $54,008 Os aie (OTe Ne >) (vh) ; 45-1008 95-1 
ZAG 
* 
ee a5 SP 85 hss }he-1 5-17Fp-1 i-1 
0 SS ree £. (uh) ; . = 
= a AO 
(5) 
Ani 5 = - (uh) 3413 si (uh) 5-13 
At a cos 0; Ad 
. (vh) 9441 COS O34) - (vh)g5-1] cos 05-] fe 
AO 
where 
= * 
Ogg = ea atie 
and 
* 
ig ~ %5"; 
where 
* S i 1 : 7 | | 
sg > 28g + Pi-ag) 7 pp {Paery 7 Big ~ Pa-ag + Pi-2;’ 


which is a second order, one-dimensional Bessel's interpola- 


tion scheme with p = . 
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Similarly 


(uh) F435 


72( (uh) 3 3+ (uh) i+13! -_ earl (uh) 44257 (uh) ct) 


= (uh) ijt (uh) i-13 ] 


and 
* ae ae | = 


-—(vh) . 


fon gaaiated 


2 G41! 


Sercginally, all the interpolated, starred quantities were 
derived using a two-dimensional linear interpolation. This 
method was found to introduce a 2d-wave which in certain 
areaS was an order of magnitude greater in amplitude than 
the zonal wave used for the experiments. 

Figure 2 shows the indexing convention used for the 
mass and wind variables. The index 1 in equations (5) and 


(6) was i if j was odd and 1+ 1 if j was even. 


Heodd -Wi4 “Mi + , 
cen . ~ + Wis . M45 
J odd “Wig Mis” . 


FIGURE 2. GRID INDEXING 
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eee eis DLE PRRENCING METHODS 


Four time differencing methods were used to evaluate the 
phase angle and amplitude errors of each method. The errors 
were evaluated by comparison to an analytic solution to the 
non-divergent barotropic vorticity equation. 

The four methods tested were the leapfrog, Euler-backward, 
leapfrog-trapezoidal, and Adams-Bashford schemes. Two tests 
were made with each method. 

The first set of tests used time increments of fifteen 
minutes for the first three methods and ten minutes for the 
Adams-Bashford scheme. 

These time steps were possible, even though the i grid 
distance at 85° N and S is only 40 km, because of a pro- 
cedure, used by Arakawa (Langlois and Kwok, 1969) to average 
the effects of high frequency inertial gravity waves in the 


zonal direction. 


The averaging technique involved a coefficient 


A. = »125(D; - 1) /D; (7) 
where 
D. = AS 
COs 05 


and D was the greatest integer value of D5. 
The uh terms in equation (3) and in the advective terms 
of equation (1) and (2) were replaced by 


(uh)? = (uh), +A .[(uh), | +(uh), | .-2(uh)..] 
Ji) sa) 3 sla AGG ioe shy 


IL) 


For l< D5 Z (8) 


Ge: 





and 


(uh)N = (uh) Nclea [ (uh) N71 4+ (uh) No} -2 (uh) N72] 
sy 13) J tae lly 4-15 sia) 


for N<D3< N+1 (9) 


Similarly h in the pressure gradient term in equation (1) 


was replaced by 


ht =h +A [hth = -2h J for 1<D.<2 (10) 
‘ies a) hei 13 aise J 
and 
fe 1N-1 Noten l Ay Nd | 
ae eee 4 in ee ee ie aun || seve NsD5<N+1 | (11) 


No averaging was done if Dj<l. 

The second set of tests were run without uSing the 
averaging technique. Thus to remain within the von Neumann 
ienear computational stability criterion (Haltiner, 1971), 
a 2.5-minute time step was used. The two sets of tests 


were run in order to determine the effects of the averaging 


technique on the solutions. 


fee LEAPFROG 


The leapfrog method is a centered time differencing 


scheme which is conditionally stable for eee <1. The finite 
Ax 
difference equation is: 
tee 
pttl = pt-l 4 gat oF (12) 
ot 


Since this method has three time levels, it has both a 


¢ 


physical and a computational mode. 
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B. EULER-BACKWARD 


The Euler-backward method is a two-step iterative 


scheme which is conditionally stable for CAt <1. The 


Ax 
difference equations are: 
, grt 
Beo= Fe + At __ 
dt 
oF” 


pret = Ft + At (193) 





oh 


Since the Euler-backward method has just two time levels 


menaas Only a physical mode. 


C. LEAPFROG-TRAPEZOIDAL 
The leapfrog-trapezoidal is another two-step iterative 


scheme which is conditionally stable for cue <2. 





Ax 
The difference equations are: 
F | art 
Pee Ferd + 2At —— 
dt 
At for* o9Ft 
perl = pe + : ce (14) 
A ie dt 


Since this method, like the leapfrog, has three time levels, 


it also has both a physical and a computational mode. 


D. ADAMS-BASHFORD 
The Adams-Bashford method used was the one examined by 
meily (1965). 


The difference equation is: 


owes Petpet | 5 oe — (15) 
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fees method has three time levels, thus it has both a com- 
putational and a physical mode. The method is unstable but 
has some desirable features. The computational mode tends 
to damp and the rate of erroneous amplification of the 


physical mode is small if At is small. 
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ieee ea COM DALL TONS 


The initial velocity and height fields for these ex- 
periments were derived from a stream function which is a 
solution to the non-divergent barotropic vorticity equa- 
tion. The stream function used was examined by Gates 
(1962) and Neamtan (1946), which is: 


y = A sin(mA-vt) ae (sinO) - B a2 sind +C eS HLT) 
(16) 


A reasonable meteorological pattern was obtained from equa- 
tion (16) by selecting 
C= 0 


A 1000 m¢* secu7l 


The constant B was related to the angular wave speed by 


Vv n(n+1)-2 222 
en ee (17) 
m n(n+1) n(n+1) 
For wave number 6 and, with n = 7 for convenience, Yves 20° 
m 


Bemg., per day 
B = 6.8905 x 107© secu7l 
The stream function then became 
y = -279.68 x 10-6sinO + 136.65 x 10-© sin(6A-vt) sind 
cos®Q@ mé¢sec~1 CES) 
Since these experiments were performed using a free sur- 
face barotropic primitive equations model which allows di- 
vergence, equation (17) was satisfied only approximately. 
Rossby (1939) has shown that the presence of divergence in 


a barotropic atmosphere will slow up the rate of wave propa- 


gation, especially for small values of wave number m. 
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The initial wind field was a non-divergent wind given 


by 
3p dL oop 
tea = = _. (3) 
0g a dO : 
3Y tl dW 
ee __ (20) 


ox a cos 0 ddA 


The initial height field was derived by solving the 


linear balance equation 





1 
V7h = _ [£V2) + Vw-vVE] (21) 
g 
where 
1 it 92 ils a , 
2 2 ae — (cos 0 =| 
a2 Lcos2 0 342 = cos O 30 30 


and 
( 1 9 1 9 ) 
V={—_— — ,- — 
a cos 0 dA a 00 


Equation (21) was solved by the following relaxation 








scheme: 
= SR. | (22) 
a a 1) 1) 
where 
9 | 
Rj 3 a (cos0.-d sin0j)h,.,,+(cos0;+d sin0j)hj5_5 
(cos0O.+ ) 
J coso: 
J 
(Cage Beane 1 cos@. (2a) 
ee (cose hyy (fy “y+ Vy VE) 
cos0. cos0 g 
(23) 


With a relaxation tolerance of .1 meters. 
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One experiment, using the leapfrog scheme for time 
differencing, was performed with the “restorative-iterative" 
initialization method developed by Winninghoff (1971). 

This method involved using the Euler-backward time dif- 
Peeencing Scheme to alternately step forward and backward 
Six times. After each iteration of equations (1), (2), 


and (3) the following restoration was added: 


(uh) ig = (AS (uh) 1g ap Kk, (uh) | 
(vh) ; = (eee) (Vh) j 5 + Ky (vh) (24) 
where the k's are functions of latitude. i and i were 


meetom Latitude 20° S to 20° N, O from 40° N and S to the 
poles, and a linear variation between 0 and .5 between 
latitude 20° and 40°. ky was 25 from 40° N and S to the 
poles, 0 between 20° S and 20° N, and a linear variation 


between 20° and 40°. 


Zul 





V. WAVE ANALYSIS METHOD 


iewcalculate the phase angles and amplitudes, a fourier 
analysis was performed at each five degrees of latitude 
around the latitude circle. The fourier series was ex- 


pressed as follows: 


F(x) = A, + >> (A, cos mx + B_, sin mx) 
m 
= Cy +2 Ch cos Cnix= cr) 
where 
B A 
m m 
lq 3 
Sin (dp) eos (6,)) 
and 
B 
= -)} m 
om = tan a 
Am 


Since the input stream function involved only wave number 
Six and a mean height, only Cor Ce, and O¢ values were ex- 


tracted from the fourier analysis. 
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VI. RESULTS 


All the experiments performed with the Arakawa averaging 
technique showed a considerable tilt backward at high lati- 
tudes in the phase propagation of the oven tis as te be 
expected since the smoothing of the gradients in the 
technigue tends to slow down the rate of propagation. 

Gates -(1959) has shown that, as the wavelength and Ax de- 
crease proportionally, the phase speed of the wave remains 
constant, and also that if the wavelength decreases and 

Ax remains constant, the phase speed will decrease relative 
to the exact value. In this model, the Arakawa averaging 
technique gives an effective Ax which is comparable to that 
at low latitudes, thus as the wavelength decreased toward 
the poles the phase speed also decreased. 

The result of this differential movement was to cause, 
Peentually, the formation of closed highs and lows at the 
higher latitudes which propagated equatorward. It is 
believed that this instability is possible due to non- 
linear effects introduced after the field ceased to be 
Marmonic in the latitudinal direction. 

The amplitudes in all the experiments showed a tendency 
to decrease at latitudes below 45° and increase above 45°. 
The mean height also tended to increase at the higher 
latitude (75° and above). These amplitude variations are 
also believed to be caused by the nonlinear effects. All 
the methods, except the Euler-backward, had small amplitude 


gravity waves propagating with about a ten hour period. 
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Table I shows the comparison of the time required for 
meee nour Forecast using each of the four metheds. It 
also gives a comparison of the initial twenty-four phase 


speed for selected latitudes. 


TABLE I 
"lise -Differencing} Time RequirediBhase Speed for Initial 24hrs 
Method 120hr Forecast |Equator | 30° | 60° | 
Leapforg a2 


Euler-Backward 57 


ifeapirog— 
Trapezoidal 58 


Adams-Bashford 46 





Note: Time in minutes 
Phase speed in degrees longitude per day 

The experiments performed without the averaging technique 
still showed a slight tendency to tilt backwards at the 
higher latitudes. This was not expected and was believed to 
have been caused by truncation errors due to special treat- 
ment near the poles using only a second-order difference 
approximation for the derivatives. This belief was based 
On the fact that the input stream function does not vary 
linearly near the poles, which caused problems earlier in 
the interpolation for (uh) * and h*. 

Table II gives a relative comparison of the time re- 


quired by each method uSing a 2.5-minute time step. 
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ALN SEE 


Time Differencing Time Required in Minutes 
Method 


Leapfrog 40 min. for a 32 hr. forecast 









Euler-Backward 70 min. for a 30 hr. forecast 






Leapfrog-Trapezoidal VOentaeetoG a,29 hr. forecast 





Adams-Bashford forecast 


A. RESULTS OF INDIVIDUAL EXPERIMENTS USING THE AVERAGING 
TECHNIQUE 


Experiment 1. This experiment was performed uSing the 
leapfrog time differencing method. Figure 3 shows the phase 
angles as a function of latitude at twelve-hour intervals 
out to thirty-six hours and Fig. 8a shows phaSe angles at 
twenty-four hour intervals out to 120 hours. The amplitudes 
for wave number six and the mean heights are shown in Fig. 
13 for selected latitudes. 

Experiment 2. The second experiment was performed using 
the leapfrog scheme and the “restorative-iterative" initial- 
ization method. This experiment was performed to see if the 
tilt of the phase lines could be reduced by letting the mass 
and wind fields "adjust" before performing the integrations. 
As can be seen in Fig. 4, the tilt was not reduced. 

Experiment 3. The Euler-backward method was used for 
this experiment. The phase curves are shown in Figs. 5 and 
8b. The mean height and wave number six amplitudes are 
Shown in Fig. 14. It should be noted that the gravity waves 
present in the other three methods are effectively damped 


Sue with this method. 
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Also the maximum variations in amplitudes, which are 
approximately equal for the other three aed - are 
slightly less since this scheme tends to damp all waves. 

' Experiment 4, This experiment was performed using the 
leapfrog-trapezoidal method. Figures 6 and 8c show the 
phase angles vs latitude curves. The amplitudes are shown 
in Fig. 15. The largest gravity wave amplitudes were ob- 
served using this method. 

Experiment 5. The last experiment using the Arakawa 
averaging scheme was performed using the Adams-Bashford 
method. Figures 7 and 8d show the phase relationships and 
Fig. 16 and amplitudes. There was very little difference 
in the results between this method and the leapfrog, except 
for the time required, see Table I, for the integrations. 
PeeeeRkEoULTS OF INDIVIDUAL EXPERIMENTS WITHOUT THE AVERAGING 

2B CHNIOUE 

Experiment 6. This experiment, like the first experi- 
ment, was performed using the leapfrog scheme. The time 
step was reduced from 15 minutes to 2.5 minutes. The phase 
angle peeeinles avemcmownela Fag, 9 for zero, twelve, and 
twenty-four Sees. 

Experiment 7. This experiment was the same as experiment 
3, except At was 2.5 minutes. The phase angle results are 
eaewn in Fig. 10. 

Experiment 8. This experiment was the pana ae experi- 
ment 4 with the exception of At, which was reduced to 2.5 


Minutes. The phase profiles are shown in Fig. ll. 
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Experiment 9. The same time differencing method was 
used as in experiment 5. The time increment was reduced 
from 10 minutes to 2.5 minutes. Figure 12 shows the phase 


relationships for this experiment. 
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3. PHASE ANGLE VS LATITUDE USING THE LEAPFROG 
SCHEME WITH THE ARAKAWA AVERAGING. 


The input height field was constant at 85° lat. The 


phase angle at that latitude is the result of trunca- 
tion‘errors in the fourier analysis. 
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SCHEME WITH THE ARAKAWA AVERAGING AND 
WINNINGHOFF'S “RESTORATIVE-ITERATIVE" 
INITIALIZATION. 
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FIGURE 5. PHASE ANGLE VS LATITUDE USING THE EULER- 
BACKWARD SCHEME WITH THE ARAKAWA AVERAGING. 


see note on Fig. 3. 
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FIGURE 6. PHASE ANGLE VS LATITUDE USING THE LEAPFROG- 
TRAPEZOIDAL SCHEME WITH THE ARAKAWA AVERAGING. 


See note on Fig. 3. 
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FIGURE 7. PHASE ANGLE VS LATITUDE USING THE ADAMS- 
BASHFORD SCHEME WITH THE ARAKAWA AVERAGING. 


See note on Fig. 3. 


32 





POunOne Oat 
OL LINO STVAXYGYLNI YNOH—-Pe LV ONIOVYAAY VMVAVEV FH 
HLIM SHWHHOS YNOT TIV YOd PAGNGLILVT SA WIONVY dSVvHd 


TWOLIOZddVad 
CYOdHSVE-SWVdV -J0U4 dV qCUYUvyMMOVd -ad ato 
opn3 thuo7 





°8 ganoid 


JOddT dVaT 


SpPnzt7eT 


a6 





Latitude 


90 





[+> 

(> 

Jo 

1 

fo] © 

| °) 

10° 20° 30° 40° 50° 60 

Longitude : 


FIGURE 9. PHASE ANGLE VS LATITUDE USING. THE LEAPFROG 
SCHEME WITHOUT AVERAGING. 


Moennoce On Fig. 3. 


This 24 hour movement appears to be wrong compared to 
the averaging case but time did not permit a re-run to 
verify this movement. 
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BPeGuURe tO. PHASE ANGLE VS LATITUDE USING THE EULER- 
BACKWARD SCHEME WITHOUT AVERAGING. 





See note on Fig. 3. 
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FIGURE ll. PHASE ANGLE VS LATITUDE USING THE LEAPFROG- 
TRAPEZOIDAL SCHEME WITHOUT AVERAGING. 


SecemOte On Fig. 3. 
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FIGURE 12. PHASE ANGLE VS LATITUDE USING THE ADAMS- 
BASHFORD SCHEME WITHOUT AVERAGING. 


See note on Fig. 3. 
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Ve CONCHUS LIONS 


The Arakawa averaging method caused some problems with 
the initial field which was represented by a spherical 
harmonic, but it is felt that with real SEA. where the 
moeieeucinal scale does not necessarily decrease with 
latitude, the method might not cause such severe problems. 
Considering the alternatives, such as a reduced time step, 
variable grid size, or variable time step, the Arawaka 
procedure 1S a simple and effective method for spherical 
prediction. The reduced time step is much too expensive 
in computer time to be practical. The abruptly changed 
grid size causes severe problems around the area of the 
change. A variable time step might prove to be acceptable 
but would involve some very complex programming. 

In the experiments performed, a second order one- 
dimensional interpolation was used since problems arose 
from using a two-dimensional linear interpolation and was 
the easiest to apply. In the real data cases, a two- 
dimensional second order interpolation scheme would prob- 
ably give better overall results. 

Overall the Euler-backward method gave the best results 
Since it effectively reduced the smh eles of the gravity 
waves, but was expensive in computer time. Considering 
time requirements and overall results, the leapfrog method 
is Still the most desirable. Some further tests with com- 
binations of the methods might produce a method which gives 
good results and is acceptable as far as time required is 
concerned, 
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